import pandas as pd
import numpy as np
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
import seaborn as sns
from sklearn import metrics
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn import preprocessing
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import plotly.graph_objs as go
import plotly.offline as offline
offline.init_notebook_mode()
Credits: based on https: // www.kaggle.com/crawford/principle-component-analysis-gene-expression/notebook
https://www.kaggle.com/crawford/principle-component-analysis-gene-expression/
Datos usados para clasificar pacientes con acute myeloid leukemia (AML) and acute lymphoblastic leukemia (ALL).
Golub et al "Molecular Classification of Cancer: Class Discovery and Class Prediction by Gene Expression Monitoring"
There are two datasets containing the initial (training, 38 samples) and independent (test, 34 samples) datasets used in the paper. These datasets contain measurements corresponding to ALL and AML samples from Bone Marrow and Peripheral Blood. Intensity values have been re-scaled such that overall intensities for each chip are equivalent.
testfile = 'genes/data_set_ALL_AML_independent.csv'
trainfile = 'genes/data_set_ALL_AML_train.csv'
labels = 'genes/genes.actual.csv'
X_train = pd.read_csv(trainfile)
X_test = pd.read_csv(testfile)
y = pd.read_csv(labels)
# 1) Remove "call" columns from training a test
train_keepers = [col for col in X_train.columns if "call" not in col]
test_keepers = [col for col in X_test.columns if "call" not in col]
X_train = X_train[train_keepers]
X_test = X_test[test_keepers]
X_train.head()
| Gene Description | Gene Accession Number | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | ... | 35 | 36 | 37 | 38 | 28 | 29 | 30 | 31 | 32 | 33 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | AFFX-BioB-5_at (endogenous control) | AFFX-BioB-5_at | -214 | -139 | -76 | -135 | -106 | -138 | -72 | -413 | ... | 7 | -213 | -25 | -72 | -4 | 15 | -318 | -32 | -124 | -135 |
| 1 | AFFX-BioB-M_at (endogenous control) | AFFX-BioB-M_at | -153 | -73 | -49 | -114 | -125 | -85 | -144 | -260 | ... | -100 | -252 | -20 | -139 | -116 | -114 | -192 | -49 | -79 | -186 |
| 2 | AFFX-BioB-3_at (endogenous control) | AFFX-BioB-3_at | -58 | -1 | -307 | 265 | -76 | 215 | 238 | 7 | ... | -57 | 136 | 124 | -1 | -125 | 2 | -95 | 49 | -37 | -70 |
| 3 | AFFX-BioC-5_at (endogenous control) | AFFX-BioC-5_at | 88 | 283 | 309 | 12 | 168 | 71 | 55 | -2 | ... | 132 | 318 | 325 | 392 | 241 | 193 | 312 | 230 | 330 | 337 |
| 4 | AFFX-BioC-3_at (endogenous control) | AFFX-BioC-3_at | -295 | -264 | -376 | -419 | -230 | -272 | -399 | -541 | ... | -377 | -209 | -396 | -324 | -191 | -51 | -139 | -367 | -188 | -407 |
5 rows × 40 columns
# 2) Transpose
X_train = X_train.T
X_test = X_test.T
X_train.head()
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 7119 | 7120 | 7121 | 7122 | 7123 | 7124 | 7125 | 7126 | 7127 | 7128 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Gene Description | AFFX-BioB-5_at (endogenous control) | AFFX-BioB-M_at (endogenous control) | AFFX-BioB-3_at (endogenous control) | AFFX-BioC-5_at (endogenous control) | AFFX-BioC-3_at (endogenous control) | AFFX-BioDn-5_at (endogenous control) | AFFX-BioDn-3_at (endogenous control) | AFFX-CreX-5_at (endogenous control) | AFFX-CreX-3_at (endogenous control) | AFFX-BioB-5_st (endogenous control) | ... | Transcription factor Stat5b (stat5b) mRNA | Breast epithelial antigen BA46 mRNA | GB DEF = Calcium/calmodulin-dependent protein ... | TUBULIN ALPHA-4 CHAIN | CYP4B1 Cytochrome P450; subfamily IVB; polypep... | PTGER3 Prostaglandin E receptor 3 (subtype EP3... | HMG2 High-mobility group (nonhistone chromosom... | RB1 Retinoblastoma 1 (including osteosarcoma) | GB DEF = Glycophorin Sta (type A) exons 3 and ... | GB DEF = mRNA (clone 1A7) |
| Gene Accession Number | AFFX-BioB-5_at | AFFX-BioB-M_at | AFFX-BioB-3_at | AFFX-BioC-5_at | AFFX-BioC-3_at | AFFX-BioDn-5_at | AFFX-BioDn-3_at | AFFX-CreX-5_at | AFFX-CreX-3_at | AFFX-BioB-5_st | ... | U48730_at | U58516_at | U73738_at | X06956_at | X16699_at | X83863_at | Z17240_at | L49218_f_at | M71243_f_at | Z78285_f_at |
| 1 | -214 | -153 | -58 | 88 | -295 | -558 | 199 | -176 | 252 | 206 | ... | 185 | 511 | -125 | 389 | -37 | 793 | 329 | 36 | 191 | -37 |
| 2 | -139 | -73 | -1 | 283 | -264 | -400 | -330 | -168 | 101 | 74 | ... | 169 | 837 | -36 | 442 | -17 | 782 | 295 | 11 | 76 | -14 |
| 3 | -76 | -49 | -307 | 309 | -376 | -650 | 33 | -367 | 206 | -215 | ... | 315 | 1199 | 33 | 168 | 52 | 1138 | 777 | 41 | 228 | -41 |
5 rows × 7129 columns
# 3) Clean up the column names for training data
X_train.columns = X_train.iloc[1]
X_train = X_train.drop(["Gene Description", "Gene Accession Number"]).apply(pd.to_numeric)
# Clean up the column names for training data
X_test.columns = X_test.iloc[1]
X_test = X_test.drop(["Gene Description", "Gene Accession Number"]).apply(pd.to_numeric)
X_train.head()
| Gene Accession Number | AFFX-BioB-5_at | AFFX-BioB-M_at | AFFX-BioB-3_at | AFFX-BioC-5_at | AFFX-BioC-3_at | AFFX-BioDn-5_at | AFFX-BioDn-3_at | AFFX-CreX-5_at | AFFX-CreX-3_at | AFFX-BioB-5_st | ... | U48730_at | U58516_at | U73738_at | X06956_at | X16699_at | X83863_at | Z17240_at | L49218_f_at | M71243_f_at | Z78285_f_at |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | -214 | -153 | -58 | 88 | -295 | -558 | 199 | -176 | 252 | 206 | ... | 185 | 511 | -125 | 389 | -37 | 793 | 329 | 36 | 191 | -37 |
| 2 | -139 | -73 | -1 | 283 | -264 | -400 | -330 | -168 | 101 | 74 | ... | 169 | 837 | -36 | 442 | -17 | 782 | 295 | 11 | 76 | -14 |
| 3 | -76 | -49 | -307 | 309 | -376 | -650 | 33 | -367 | 206 | -215 | ... | 315 | 1199 | 33 | 168 | 52 | 1138 | 777 | 41 | 228 | -41 |
| 4 | -135 | -114 | 265 | 12 | -419 | -585 | 158 | -253 | 49 | 31 | ... | 240 | 835 | 218 | 174 | -110 | 627 | 170 | -50 | 126 | -91 |
| 5 | -106 | -125 | -76 | 168 | -230 | -284 | 4 | -122 | 70 | 252 | ... | 156 | 649 | 57 | 504 | -26 | 250 | 314 | 14 | 56 | -25 |
5 rows × 7129 columns
# 4) Split into train and test
X_train = X_train.reset_index(drop=True)
y_train = y[y.patient <= 38].reset_index(drop=True)
# Subet the rest for testing
X_test = X_test.reset_index(drop=True)
y_test = y[y.patient > 38].reset_index(drop=True)
# Obtain numerical labels in train and test
le = preprocessing.LabelEncoder().fit(y_train["cancer"])
# To numerical
y_train_num = le.transform(y_train["cancer"])
y_test_num = le.transform(y_test["cancer"])
Realiza un análisis exploratorio de los datos (correlaciones entre sí y con las clases, distribuciones,...). Usa las técnicas y gráficos que te parezcan más representativos.
Podemos comenzar mirando qué clases tenemos y la distribución que tienen las mismas
df = pd.DataFrame(y)
print(df["cancer"].unique())
sns.histplot(df["cancer"])
['ALL' 'AML']
<AxesSubplot:xlabel='cancer', ylabel='Count'>
Como podemos observar, hay más pacientes padecen ALL que pacientes que padecen AML en nuestro dataset (la proporción orbita entorno a $\frac{2}{3}:\frac{1}{3}$ respectivamente). A continuación, de cara a estudiar si existen correlaciones entre las distintas variables podemos representar un heatmap de estas correlaciones. Puesto que tenemos más de $7000$ variables, presentamos únicamente las $50$ primeras.
n_features = 50
data = X_train[X_train.columns[:n_features]].copy()
data['cancer'] = y['cancer']
# Compute the correlation matrix and triangle mask
cor_matrix = data.corr()
mask = np.triu(np.ones_like(cor_matrix, dtype=bool))
# Generate a custom colormap
cmap = sns.diverging_palette(150, 275, s=80, l=55, n=9)
# Draw the heatmap with the mask and correct aspect ratio
plt.subplots(figsize=(11, 9))
_ = sns.heatmap(cor_matrix, mask=mask, cmap=cmap, vmax=.3, center=0,
square=True, linewidths=.5, cbar_kws={"shrink": .5})
Se puede apreciar en el gráfico anterior que encontramos correlaciones de hasta $0.3$, y correlaciones negativas de casi $-0.6$. Debido a la alta dimensionalidad de los datos es díficil hacerse a la idea de las correlaciones realmente presentes en el dataset, aunque esta representación nos indica que son no triviales.
A continuación estudiamos cómo se distribuyen las dos clases en cada variable individualmente. Para ello pintamos boxplots de las primeras $50$ variables utilizando seaborn.
data = X_train.copy()
data['cancer'] = y_train['cancer']
fig = plt.figure(figsize = (20, 50))
for i, col_name in enumerate(X_train.columns[:n_features]):
plt.subplot(10, 5, i+1)
ax = sns.boxplot(x='cancer', y=col_name, data=data)
fig.tight_layout()
Podemos apreciar como las clases no son claramente separables en ninguna variable, siendo especialmente relevantes las últimas tres en las que estamos relativamente cercanos a este hecho, pero aún así una cantidad no trivial de datos (entre el $25$ y el $50\%$) se solapan a simple vista. Era de esperar que ante un problema tan complejo no pudiesemos distinguir las clases utilizando unicamente una variable.
De nuevo, este análisis no es definitivo pues hemos estudiando solo las primeras $50$ variables. Pasemos a estudiar cómo afectan los algoritmos de PCA y LDA a nuestro dataset.
The analysis reveals that 21 principle components are needed to account for 80% of the variance. PC 1-3 add up to about ~33% and the rest is a slow burn where each component after PC8 contributes between 1-2% of the variance up until PC38 which is essentially zero. 1% is a decent amonut of variance and so the number of important PCs is up for interpretation.
# 5) Scale data
# (1) YOUR CODE HERE: Use the StandardScaler (separately for train and test sets)
X_train = StandardScaler().fit_transform(X_train)
X_test = StandardScaler().fit_transform(X_test)
# 6) PCA Analysis and projection
components = 21
# YOUR CODE HERE:
# (2) Use PCA with this number of components on train set, with Y the result of the procedure
pca = PCA(n_components = components)
X_train_pca = pca.fit_transform(X_train)
# (3) Retrieve the explained variance ratio, and compute its accumulative sum
# save those values in variables var_exp and cum_var_exp
var_exp = pca.explained_variance_ratio_
cum_var_exp = np.cumsum(var_exp)
print(var_exp)
print(cum_var_exp)
[0.14987793 0.11977811 0.06600568 0.04884912 0.04632413 0.03721924 0.03490808 0.03289659 0.0298473 0.02644449 0.0250942 0.02356641 0.0220405 0.02084695 0.01935673 0.01892268 0.01840658 0.01708862 0.01700553 0.01640067 0.01529196] [0.14987793 0.26965604 0.33566172 0.38451084 0.43083497 0.46805422 0.5029623 0.53585888 0.56570619 0.59215068 0.61724488 0.64081129 0.66285179 0.68369875 0.70305548 0.72197815 0.74038473 0.75747335 0.77447888 0.79087955 0.80617152]
Pregunta (1): ¿Qué pauta puede observarse en los valores de var_exp? ¿Cuál es la interpretación relativa de esos valores?
# Plot the explained variance using var_exp and cum_var_exp
x = ["PC%s" %i for i in range(1,components)]
trace1 = go.Bar(
x=x,
y=list(var_exp),
name="Explained Variance")
trace2 = go.Scatter(
x=x,
y=cum_var_exp,
name="Cumulative Variance")
layout = go.Layout(
title='Explained variance',
xaxis=dict(title='Principle Components', tickmode='linear'))
data = [trace1, trace2]
fig = go.Figure(data=data, layout=layout)
offline.iplot(fig)
Se puede observar que a medida que avanzamos en el número de componente, el porcentaje de la varianza que explica esta componente disminuye. Esto tiene sentido, pues son las primeras componentes de las componentes principales las que son capaces de explicar la mayor variabilidad de los datos.
La variabilidad acumulada explicada, aún así, sigue aumentando cuando vamos avanzando en el número de componente, pues aunque la cantidad de varianza explicada sea mucho menor que en las componentes iniciales, sigue sin ser cero (a partir de la décima componente, todas explican menos de un $3\%$ de la varianza total, pero todas se mantienen por encima del $1\%$.
Recordamos también que el máximo número de componentes principales que podríamos tomar es igual al mínimo entre el número de características de un elemento del conjunto de datos de entrenamiento y el número de datos que tenemos en el conjunto de train, en este caso $38$, y conseguiríamos explicar así toda la varianza del conjunto.
Las tres primeras componentes solamente explica un $33\%$ de la varianza, pero vamos a obtener una representación visual de la proyección a estas tres primeras componentes.
# Project first three components
n_components = 3
pca_3 = PCA(n_components = n_components)
X_train_pca = pca_3.fit_transform(X_train)
traces = []
for name in ['ALL', 'AML']:
trace = go.Scatter3d(
x=X_train_pca[y_train.cancer == name, 0],
y=X_train_pca[y_train.cancer == name, 1],
z=X_train_pca[y_train.cancer == name, 2],
mode='markers',
name=name,
marker=go.Marker(size=10, line=go.Line(width=1), opacity=1))
traces.append(trace)
layout = go.Layout(
xaxis=dict(title='PC1'),
yaxis=dict(title='PC2'),
title="Projection of First Three Principle Components"
)
data = traces
fig = go.Figure(data=data, layout=layout)
offline.iplot(fig)
c:\users\jose\appdata\local\programs\python\python37\lib\site-packages\plotly\graph_objs\_deprecations.py:385: DeprecationWarning: plotly.graph_objs.Line is deprecated. Please replace it with one of the following more specific types - plotly.graph_objs.scatter.Line - plotly.graph_objs.layout.shape.Line - etc. c:\users\jose\appdata\local\programs\python\python37\lib\site-packages\plotly\graph_objs\_deprecations.py:441: DeprecationWarning: plotly.graph_objs.Marker is deprecated. Please replace it with one of the following more specific types - plotly.graph_objs.scatter.Marker - plotly.graph_objs.histogram.selected.Marker - etc.
Pregunta(2): Modificando la perspectiva de la figura con el ratón, ¿qué observas en cuanto a la separabilidad de las clases? Adjunta una imagen que apoye tus conclusiones.
Realizamos dos perspectivas de los datos que nos pueden dar una idea intuitiva de lo que las tres primeras componentes principales nos dicen sobre los daots.
Si giramos la figura del siguiente modo y trazamos una línea negra que representa un posible hiperplano de separación obtenemos:

Vemos como parece que los puntos de la clase roja AML parecen estar más a la izquierda y por encima que los de la clase azul ALL.
Girando más la perspectiva, de forma que la clase roja quede a la derecha, vuelve a parecer que hay un grado alto de separabilidad en nuestros datos.

Realizar un análisis similar usando LDA, usando en este caso la información sobre el tipo de cancer de cada paciente.
# Define LDA model
lda = LDA(n_components = 1)
# Fit and transform Train data
X_train_lda = lda.fit_transform(X_train, y_train["cancer"])
def plot_lda_projection(X, y, title, cluster_centers=None):
labels = ['ALL', 'AML']
markers = ['s','o']
colors = ["orange", "blue"]
ax = plt.subplot()
for idx, (label, marker, color) in enumerate(zip(labels, markers, colors)):
if cluster_centers is None:
cc = np.mean(X[y == idx])
else:
cc = cluster_centers[idx]
plt.scatter(x=X[y== idx], y=np.zeros(len(X[y==idx])),
marker=marker,
color=color,
alpha=0.5,
label=label)
plt.plot(cc, 0.0,
marker = "x",
color = color,
markersize = 14,
label = "Mean of {}".format(labels[idx]))
plt.xlabel('Projection1')
plt.legend(loc='upper right', fancybox=True)
plt.title(title)
ax.set_yticks([0])
for pos in ['top','right','bottom','left']:
ax.spines[pos].set_visible(False)
Dibujamos primero un gráfico para observar el resultado de la proyección a $1$ dimensión usando LDA.
plot_lda_projection(X_train_lda, y_train_num, "LDA Projection Train")
Como podemos comprobar, tras proyectar los datos usando linear discriminant analysis, obtenemos que los elementos de las clases son prácticamente linealmente separables, considerando que puede haber datos que sean ruido (algunos elementos de la clase $0$ se clasificarían como los de la clase $1$).
Utiliza k-means clustering con los datos originales y con los datos proyectados con PCA y LDA. ¿Qué observas?
Recordamos que, hasta este punto, tenemos:
X_train el conjunto de datos de entrenamiento estandarizadosX_train_pca el conjunto de entrenamiento al que se le ha aplicado PCA con $n=3$ componentesX_lda_sklearn el conjunto de entrenamiento al que se le ha aplicado LDA con $n = 1$ componenteUsualmente, Kmeans es un algoritmo de clustering, no de clasificación. Sin embargo, en este caso podemos tratar de enfocar el problema de la siguiente manera:
KMeans usando $k = 2$, que es el número de clases.Pasamos a realizar el análisis usando KMeans.
Recordamos también que este modelo tiene como función de pérdida la inertia o suma de las distancias entre los elementos del mismo cluster, que puede expresarse como
$$ \sum_{i=0}^{n}\min_{\mu_j \in C}||x_i - \mu_j||^2. $$Hay que tener en cuenta que es una métrica no normalizada, que nos indica cómo de coherentes son los clusters que se han creado.
Tenemos que buscar también métrica de evaluación. Claramente podemos usar esta inertia como una de estas métricas, pero nos interesa buscar alguna más. En este caso, sabemos las clases originales que tienen las muestras, por lo que podemos buscar una métrica de external-cluster validation. Nos interesa una que evalúe las muestras dos a dos, así que usaremos el Rand Index, en concreto, el adjusted Rand Index. En general, esta es una medida de cuánto se parecen dos conjuntos de clusters sobre un mismo dataset. En nuestro caso, compararemos el clustering generado por KMeans con el 'clustering' generado al asociar a cada elemento, su clase. Si estos dos conjuntos de clusters son parecidos, habremos obtenido una buena distribución.
Formulemos primero el modelo no ajustado y luego lo ajustaremos. Sea $C$ el conjunto de clustering heredado de las clases del problema y $K$ el clustering definido sobre los datos. Definimos
Si $C_2^{n_\text{samples}}$ es el número total de posibles parejas del conjunto de datos, el unadjusted Rand index se define entonces como:
$$ RI = \frac{a + b}{C_2^{n_\text{samples}}}. $$Sin embargo, para garantizar que asignar etiquetas aleatorias a los puntos nos dé un $RI$ cercano a cero, tenemos que ajustar este coeficiente definiendo el adjusted Rand index como sigue:
$$ ARI = \frac{RI - \mathbb E[RI]}{max(RI) - \mathbb E [RI]} $$Este coeficiente es simétrico y está en el intervalo $[-1,1]$, siendo mejor el clústering cuando este coeficiente se acerca a $1$ y peor cuando se aleja de este.
Recordamos que en estas métricas, lo que queremos es minimizar la inertia, y maximizar el ARI.
Creamos primero una función que ajuste $k = 2$ medias usando el conjunto de datos y nos devuelva las etiquetas predichas para ambos conjuntos de train y test.
def kmeans_predict(X_train, X_test):
# Initialize and fit kmeans
kmeans = KMeans(n_clusters = 2, random_state = 0)
kmeans.fit(X_train)
# Predict
y_train_pred = kmeans.predict(X_train)
y_test_pred = kmeans.predict(X_test)
return y_train_pred, y_test_pred, kmeans.inertia_, kmeans
def plot_results(train_pred, test_pred, train_true, test_true, inertia, name):
print(name)
print("\t + Inertia : {}".format(inertia))
print("\t + Adjusted Rand Index")
print("\t\t - Train : {}".format(metrics.adjusted_rand_score(train_true, train_pred)))
print("\t\t - Test : {}".format(metrics.adjusted_rand_score(test_true, test_pred)))
# tro stands for train original, teo stands for test original
tro_pred, teo_pred, ori_inertia, _ = kmeans_predict(X_train, X_test)
# Results
plot_results(tro_pred, teo_pred, y_train_num, y_test_num, ori_inertia, name ="Original data")
Original data + Inertia : 241408.8610451321 + Adjusted Rand Index - Train : 0.15820603069777567 - Test : -0.02674123029994916
ANÁLISIS
Ahora, queremos aplicar primero PCA para después hacer $k=2$ clústers y ver cómo de bien se agrupan nuestros datos. En principio, hemos dicho que con $n=21$ variables explicamos una gran parte de la varianza de nuestro conjunto de datos, por lo que comenzaremos haciendo PCA y quedándonos con las primeras $21$ variables para ver el resultado.
# Use pca to obtain test reduced
X_train_pca = pca.transform(X_train)
X_test_pca = pca.transform(X_test)
# trp stands for train pca, tep stands for test pca
trp_pred, tep_pred, pca_inertia, _ = kmeans_predict(X_train_pca, X_test_pca)
# Results
plot_results(trp_pred, tep_pred, y_train_num, y_test_num, pca_inertia, name ="PCA data")
PCA data + Inertia : 190705.89237483937 + Adjusted Rand Index - Train : 0.20112339486030617 - Test : -0.030352138186116567
Vemos que usando $n=21$ componentes obtenemos un ligero descenso en la métrica de inertia, lo cual es positivo. Sin embargo, la métrica ARI desciende en ambos conjuntos. Como la métrica ARI es bastante relevante pues, de algún modo, nos mide el accuracy que tienen nuestros clústers, deberíamos comenzar quedándonos con realizar los clústers con el conjunto de datos original.
A continuación, realizaremos un análisis más exhaustivo de cómo varían las métricas según el número de componentes que tomemos para PCA.
# Experiment to compute ARI in train in test and inertia for n_components in PCA
max_components = min(X_train.shape[0], X_train.shape[1])
aris_train = []
aris_test = []
inertias = []
for i in range(1,max_components+1):
pca_N = PCA(n_components = i )
X_train_pca_i = pca.fit_transform(X_train)
X_test_pca_i = pca.transform(X_test)
tr_pred_i, te_pred_i, inertia_i, _ = kmeans_predict(X_train_pca_i, X_test_pca_i)
inertias.append(inertia_i)
aris_train.append(metrics.adjusted_rand_score(y_train_num, tr_pred_i))
aris_test.append(metrics.adjusted_rand_score(y_test_num, te_pred_i))
xs = np.arange(1, max_components+1, 1)
fig, ax = plt.subplots(1,3, figsize =(20,8))
ax[0].plot(xs, inertias)
ax[0].set_title("Inertias")
ax[0].axvline(x=3, linestyle= "--" ,color = "green")
ax[0].axvline(x=19, linestyle= "--" ,color = "red")
ax[1].plot(xs, aris_train)
ax[1].axvline(x=3, linestyle= "--" ,color = "green")
ax[1].axvline(x=19, linestyle= "--" ,color = "red")
ax[1].set_title("ARI in train")
ax[2].plot(xs, aris_test)
ax[2].set_title("ARI in test")
ax[2].axvline(x=3, linestyle= "--" ,color = "green")
ax[2].axvline(x=19, linestyle= "--" ,color = "red")
plt.show()
Se puede observar cómo para cualquier número de componentes, el ARI en el conjunto de entrenamiento oscila en el interalo $[0.07, 0.20]$, alcanzando su máximo en $n=19$ componentes principales. En este punto se alcanza también un máximo en el gráfico de las inertias, que recordamos que queríamos minimizar, por lo que no consideramos $n=19$ componentes principales como un buen número de componentes para realizar nuestro análisis.
Sin embargo, observando la línea verde, vemos como con $n=3$ componentes principales obtenemos un mínimo local cercano al mínimo global de la función de inertias, obteniendo también un máximo global (sin considerar el punto $n=19$) de la función ARI en el conjunto de train, por lo que posiblemente quedarnos con $n=3$ componentes principales sería más adecuado.
Para hacernos una idea geométrica del resultado tras aplicar KMeans y PCA, dibujémos las tres primeras componentes principales de los puntos, así como los centros de los clusters y las debidas clasificaciones.
# Create data for n=3 components
X_pca_3 = pca_3.transform(X_train)
X_test_pca_3 = pca_3.transform(X_test)
y_train_pred, y_test_pred, _ , model = kmeans_predict(X_pca_3, X_test_pca_3)
xs = [el[0] for el in X_pca_3]
ys = [el[1] for el in X_pca_3]
zs = [el[2] for el in X_pca_3]
cs = model.cluster_centers_
print(model.cluster_centers_)
train_correct = y_train_pred == y_train_num
test_correct = y_test_pred == y_test_num
c0_train = np.where(y_train_num == 0)[0]
c1_train = np.where(y_train_num == 1)[0]
[[-25.25664574 6.60730097 -2.38232805] [ 28.06293972 -7.34144552 2.64703116]]
markers = {'0': 'o', '1':'s'}
colors = [{'True' : 'blue', 'False' : 'orange'},{'True': "red", 'False' : "green"}]
fig = plt.figure(figsize=(14,8))
ax = plt.axes(projection='3d')
for i in range(X_pca_3.shape[0]):
marker = markers[str(y_train_num[i])]
color = colors[y_train_num[i]][str(train_correct[i])]
ax.scatter(X_pca_3[i][0], X_pca_3[i][1], X_pca_3[i][2], marker = marker, color = color,s = 50)
ax.scatter(cs[:,0], cs[:,1], cs[:,2], marker = "*", s = 1000, label = "Cluster centers" , color = "black")
blue_patch = mpatches.Patch(color='blue', label='Class 0, Well Classified')
orange_patch = mpatches.Patch(color='orange', label='Class 0, Wrongly Classified')
red_patch = mpatches.Patch(color='red', label='Class 1, Well Classified')
green_patch = mpatches.Patch(color='green', label='Class 1, Wrongly Classified')
green_patch = mpatches.Patch(color='green', label='Class 1, Wrongly Classified')
star_patch = mpatches.Patch(color='black', label='Cluster centers')
ax.legend(handles=[blue_patch, orange_patch, red_patch, green_patch, star_patch])
ax.grid(False)
#ax.set_axis_off()
Se puede observar como tenemos puntos de la clase cero (circulares) que quedan a ambos lados del centro del clúster más cercano, por lo que todos ellos están mal clasificados. Sin embargo, los de la clase 1 quedan también en su mayoría más cercanos a este mismo centro, por lo que quedan todos bien clasificados salvo 1 (en verde), que está más alejado.
Realizamos un proceso análogo de KMeans sobre los resultados obtenidos con LDA.
# Define LDA model
lda = LDA(n_components = 1)
# Fit and transform Train data and test data
X_train_lda = lda.fit_transform(X_train, y_train["cancer"])
X_test_lda = lda.transform(X_test)
# trp stands for train pca, tep stands for test pca
trp_pred, tep_pred, lda_inertia, model = kmeans_predict(X_train_lda, X_test_lda)
plot_results(trp_pred, tep_pred, y_train_num, y_test_num, lda_inertia, name ="LDA data")
LDA data + Inertia : 21.357177393370193 + Adjusted Rand Index - Train : 0.5063013138909815 - Test : 0.05815018315018315
Podemos apreciar un enorme descenso en el valor de inertia, y un aumento notable en ambos valores de ARI. Esto revela una notable mejoría en el clustering realizado respecto a utilizar PCA. De nuevo, realizamos una representación gráfica de los resultados para obtener una idea intuitiva de las separaciones obtenidas, añadiendo una línea de separación negra en el punto donde partimos nuestra clasificación entre ALL y AML.
cluster_centers = np.array(model.cluster_centers_).flatten()
division_point = np.mean(cluster_centers)
plot_lda_projection(X_train_lda, y_train_num, "LDA Projection in Train", cluster_centers)
plt.vlines(division_point, ymin=-.5, ymax=.5, color='black')
plt.show()
plot_lda_projection(X_test_lda, y_test_num, "LDA Projection in Test", cluster_centers)
plt.vlines(division_point, ymin=-.5, ymax=.5, color='black')
plt.show()
Podemos apreciar una notable mejoría en la separación realizadas en train, donde la mayoría de elementos están asociados al cluster con su clase. Es por ello que obtenemos un valor tan elevado en el ARI de training. Sin embargo, estos resultados no se traducen en una clara mejoría en test: el modelo no generelizada bien.
Como conclusión final podemos asegurar que ni PCA ni LDA nos permiten generalizar bien a partir del dataset de training, vista la enorme diferencia del ARI entre training y test en ambos casos. Esto puede deberse al bajo número de muestras en nuestro conjunto de datos. Sin embargo, utilizar LDA ha garantizado resultados notablemente mejores tanto en training como test.